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Abstract 

In order to solve the magnetohydrodynamics (MHD) equations with a H(div)-conforming 
element, a novel approach is proposed to ensure the exact divergence-free condition on the 
magnetic field. The idea is to add on each element an extra interior bubble function from 
higher order hierarchical 'H(div)-conforming basis. Four such hierarchical bases for the 7^(div)- 
conforming quadrilateral, triangular, hexahedral and tetrahedral elements are either proposed 
(in the case of tetrahedral) or reviewed. Numerical results have been presented to show the linear 
independence of the basis functions for the two simplicial elements. Good matrix conditioning 
has been confirmed numerically up to the fourth order for the triangular clement and up to the 
third order for the tetrahedral element. 
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1 Introduction 

The magnetohydrodynamics (MHD) equations describe the dynamics of a charged system under 
the interaction with a magnetic field and the conservation of the mass, momentum and energy for 
the plasma system. Such a dynamics is considered constrained as the magnetic field of the system 
is evolved with the constraint of zero divergence, namely, V-B = 0. Numerical modeling of plasmas 
has shown that the observance of the zero divergence of the magnetic field plays an important role 
in reproducing the correct physics in the plasma fluid [l]. Various numerical techniques have been 
devised to ensure the computed magnetic field to maintain divergence- free [2] . In the original work 
of [l] a projection approach was used to correct the magnetic field to have a zero divergence. 

A more natural way to satisfy this constraint is through a class of the so-called constrained 
transport (CT) numerical methods based on the ideas in [3]. As noted in [4], a piecewise ■H(div) 
vector field on a finite element triangulation of a spatial domain can be a global 'H(div) field if 
and only if the normal components on the interface of adjacent elements are continuous. Thus, in 
most of the CT algorithms for the MHD, the surface averaged magnetic fiux over the surfaces of a 
3-D element is used to represent the magnetic field while the volume averaged conserved quantities 
(mass, momentum, and energy) are used. 

In the two seminal papers O [6], Nedelec proposed to use quantities (moments of tangential 
components of vector fields) on edges and faces to define the finite dimensional space in H(div) 
and 7^ (curl). The specific construction of the basis functions in both spaces, specifically in ?^(div), 
can be done in various ways such as the hierarchical type of basis proposed in [7] for both ■H(div) 
and ?^(curl). Unfortunately, the proposed hierarchical basis for ^^(div) in [7] turns out to be 
erroneous as it can be easily checked that for quadratic polynomial approximation the proposed 
edge-based basis functions happens to be linearly dependent. 
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In this paper, we will first present a new hierarchical basis functions in ^^(div) in tetrahedral, 
for completeness, together with a review of the hierarchical basis of T-L{div) for rectangles in 2-D 
and for cube in 3-D. An important common feature of the hierarchical basis functions is the fact 
that for order p > 3 in the case of simplexes, the basis functions will include interior bubble basis 
functions which have zero normal components on the whole boundary of the element boundaries. 
Therefore, a simple way to enforce the divergence-free on each element can be easily accomplished 
by adding one single {p + l)-th order (any q-th order, q > max(m,p + = 4 for tetrahedral 

element, m = 3 for triangular element, and m = 2 for quadrilateral or hexahedral element) interior 
bubble function to a p-th order ?^(div) basis. This extra bubble basis will be able to satisfy the 
local divergence-free condition. Due to the fact of normal continuity of the p-th order basis across 
the element interface and zero normal component of the added-in p + 1-th order interior bubble 
functions, the augmented function space satisfies divergence- free condition globally. 

The rest of the paper is organized as follows. The constructions of the locally divergence- 
free bases are given in Section 2-5 for rectangular and triangular elements in 2-D, and cubic and 
tetrahedral elements in 3-D. The divergence- free condition is discussed in Section 6. Numerical 
results on the matrix conditioning are given in Section 7. Concluding remarks are given in Section 
8. 

2 Basis functions for the quadrilateral element 

In [8] Zaglmayr gave a hierarchical basis for quadrilateral ?^(div)-conforming element. In this 
section we summarize the result in [8]. 

The basis functions are constructed on the reference element - a unit square Q := [0, 1]^ [8] 
with vertexes of Vi(0, 0), ^2(1,0), V3(l,l) and V4(0, 1). The coordinate system for the reference 
element is in terms of the variables (^, 77). A bilinear function Aj which is associated with a specific 
vertex Vi has been utilized for the construction. The multiplicative factors of the bilinear function 
Aj have been used to form the linear function cjj, viz. 
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The bilinear function has this favorable property 

Ailv,- = Sij, (2) 

where 6ij is the Kronecker delta. The edge e := [Vi, Vj] which points from vertex Vi to vertex Vj is 
parameterized by 

Ce:=<rj-ai£[-l,l]. (3) 
For convenience of basis construction, the linear edge-extension parameter is also defined, viz. 

Ae := Ai + Aj G [0,1], (4) 
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which is one on edge e := [V^i,V^] and zero on the opposite edge. Note that the unit tangential 
vector Te and the outward unit normal vector can be deduced as 

fe = ^VCe, ne = VA^. (5) 

2.1 Edge-based functions 

These functions are further grouped into two categories: the lowest-order and higher-order func- 
tions. 

Lowest-order functions : 

These functions are associated with the four edges. By construction each function is perpendic- 
ular to the associated edge and has unit normal component on the associated edge. Furthermore 
the divergence of each function is unit. The shape function is given by 

V'eT° = ^Ae, (V X CeJ , i = {1, 2, 3, 4}, (6) 

which has the property 

f,,-Ve1^0=0, ne,-<^«|e, = 1, V • ^.^^^ = 1. (7) 

Higher-order functions: 

The function is taken to be a curl field of a scalar function in order to be free of divergence. 
The basis function is 

= V X (CeJ) , i = {1,2,3,4}, 0<j<p-l, (8) 

where the function is the so-called integrated Legendre polynomial of degree n [8]. The 

obvious property is 

V-4+i=0, i = {l,2,3,4}, 0<j<p-l. (9) 

2.2 Interior functions 

The interior functions are classified into three categories. 
Type 1: (curl field) 

The shape functions are given as 

ipf/ = V X (Li+2(2^ - l)Lj+2{2rj - 1)) , < i, j < p - 1. (10) 
These functions are divergence free, viz. 

v-vg'=o, o<ij<p-i. (11) 
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Type 2 : 

The formula of these functions is given as 

= Li+2m - l)^;+2(2r? - l)k + ^[+2(2^ - l)L,+2(27? - l)j%, 0<i,j<p-l. (12) 
On the boundary of the reference element the normal component of these functions vanishes, viz. 

ne■^pf\e = 0, 0<i,j<p-l, (13) 
which justifies to be one type of interior functions. 
Type 3 : 

The formula of these functions is 

= L,+2(2e - = Li+2{27] - 1%,, 0<i<p-l. (14) 

Again on the boundary of the reference element the normal component of these functions vanishes, 
viz. 

ne-Vf'|e = 0, ne-^fhe = 0, < i < p - 1. (15) 

3 Basis functions for the hexahedral element 

The reference element is defined for a unit cube T-L := [0, 1]'^ in [8]. The vertexes of the cube are 
Fi(0,0,0), 1^2(1,0,0), 1^3(1,1,0), ^4(0,1,0), ^5(0,0,1), ^6(1,0,1), 1^7(1,1,1), and ^8(0,1,1). The 
basis functions are expressed in terms of the trilinear function Xj, which is one at the vertex Vj 
and zero at all other vertexes. Along with the linear function aj and in terms of the coordinates 
variables rj, Q they are given as 



(1 


-0(1 


-77)(1 


-C), 


0-1 


A2 


:=e(l 


-^)(1 


-C), 


0-2 




A3: 


= er?(i 


-C), 


0-3 


A4 


:= (1- 




-0, 


0-4 


As 


:=(1- 


-0(1- 


-v)C, 


0-5 




Ae: 




-vK, 


0-6 






A7 := 




aj 




As : 


= (1- 




0-8 



The edge e := [Vi, Vj] which points from vertex Vi 



= (l-0 + (l-r?) + (l-C), 

= e + (i-»?) + (i-C), 

= e + ^ + (i-c), 
= (i-e) + ^? + (i-c). 

= (i-0 + (i-^?) + C, 

= e + (i-??) + c, 

= ^ + v + C, 

= {l-0 + r] + C- (16) 
;o vertex Vj is parameterized by 



He ■= CFj - e [-1, 1]. (17) 

The tangential vector associated with the edge e is given by Tg = ^V/ie- The edge extension 
parameter Ag := Aj + \j E [0,1] is one on the edge e and zero on all other edges that are parallel to 



4 



the edge e. The face / = [T^j, Vk, Vi] where the vertexes Vi and Vk are not connected by an edge 
can be parameterized by 

{Cf,Vf) ■■= iai-aj,ai-ae) G [-1,1] x [-1,1]. (18) 

The hnear face extension parameter A j = Aj + Aj + + is equal to one on the face / and zero 
on the opposite face. The outward unit normal vector of face / can be obtained by rij = VAj. 

3.1 Face-based functions 

In this subsection we record the results in [8]. We have also fixed one error in [8]. These functions 
are associated with the six faces whose formulas are classified into two groups. 

Lowest-order Raviart-Thomas functions : 

^PY''=Xfnf, i = l,2,---,Q. (19) 

Higher-order functions: (divergence-free) 

These functions are constructed as curl fields of certain vectors in order to be divergence-free. 
The formulas are given as [8] 

V'/^- = V X (A^ (Lj+2(??/)VLi+2(e/) - L*+2(C/)VLj+2(7?/))) , 0<i,j<p-l, /c = l,2,--- ,6. 

(20) 

= V X (A/Li+2(e/)Vr/^), 0<i<p-l, k = l,2,---,6. (21) 
V-/^ =Vx (A/L,+2(??/)Ve/), 0<j<p-l, k = l,2,---,6. (22) 

3.2 Interior functions 

The interior functions are further classified into three categories. The triplet (Cij^iCs) •= (2^ ~ 
1, 2r] — 1, 2C — 1) is used in the formulas. The function ^„(») is the classical un-normalized Legendre 
polynomial of degree n. While we here record the results in [8], we have implemented the correction 
of a number of mistakes in [8] as well. 

Type 1: (divergence- free) 

These functions are taken to be the curl fields of certain vector functions that are associated 
with the ■H(curl)-conforming element. 

4£i+i(ei)Vi(%)^fc+2(C3)kC' < i,j,k<p-l. 

- 2L,+2(??2)4+l(C3)jr„ 0<j,k<p-l. 



v^,';;^, = 4L,+2(ei)^i+i(?y2)4+i(C3)i« - 

V'S = 4^.+i(ei)^,+2(r?2)4+i(C3)j% - 
V^J^I =2£,+i(r?2)Lfc+2(C3)kc- 
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= 2Li+2(ei)4+i(C3)i5 - 2£,+i(ei)Lfc+2(C3)kc. < i, A: < p - 1. 
V^J^/ = 2L,+2(ei)^i+i(??2)i? - 2^,+i(6)L,+2(r?2)j%, < J < p - 1. (23) 

Type 2 : 

These functions are linear combinations of certain components in the above type. 

'^i,j,k = ^i+2(6)^i+l(f?2)4+l(C3)ic +4+l(6)^i+2(r?2)4+l(C3)jr,, < i,j,k < p - I. 
'^jl = ^j+2(??2)4+l(C3)j») + ^j+l(r/2)ifc+2(C3)kc, <j,k <p-l. 

V^jf = £i+i(6)Lfc+2(C3)kc + L,+2(ei)4+i(C3)i5, 0<i,k<p-l. 

= Li+2{Cl)ij+l{V2)k + ^i+l(ei)ii+2(%)jr„ 0<i,j<p-l. (24) 

Type 3 : 

These functions are taken as certain components in Type 2. 

V^f ^ = L,+2(ei)i5, 0<i<p-l. 

= Lk+2iC3)h^ 0<k<p-l. (25) 

Two remarks are in place. 

• All the interior basis functions are linearly independent, which can be verified easily. 

• The normal traces of these interior functions vanish on the boundary dJi of the reference 
hexahedral element. This is due to the fact that on a certain face / either one of the standard 
unit vectors is perpendicular to the normal vector of the face Uf or to the fact that the 
integrated Legendre polynomials are evaluated at 1 or —1, which is zero. 

4 Basis functions for the triangular element 

The result on the basis construction has been reported in [9]. For the completeness of the current 
study, we record the basis functions in this section. 

Any point in the 2-simplex is uniquely located in terms of the local coordinate system {£,,1]). 
The vertexes are numbered as vo(0, 0), vi(l, 0), V2(0, 1). The barycentric coordinates are given as 

Ao:=l-e-r?, Ai:=e, X2:=V- (26) 

The directed tangent on a generic edge = [ji,j2] is similarly defined as in (I40p for the three- 
dimensional case. In the same manner the edge is also parameterized as in (j4ip . A generic edge 
can be uniquely identified with 

ej ■■= [ji, j2], ii = {0, 1}, ii < i2 < 2, j= ji + j2. (27) 
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The two-dimensional vectorial curl operator of a scalar quantity, which is used in our construction, 
needs a proper definition. We use the one given in the book [11], viz. 



curl(ti) := V X n : = 



du du 



(28) 



Based upon the newly created shape functions for the three-dimensional '?^(div)-conforming tetra- 
hedral elements and using the technique of dimension reduction we construct the basis for the 
7^(div)-conforming triangular elements in two dimensions. However, it is easy to see that the two 
groups for the face functions cannot be appropriately modified for our purpose. Instead we borrow 
the idea of Zaglmayr in the dissertation [8], viz., we combine the edge-based shape functions in [8] 
with our newly constructed edge-based and bubble interior functions. In [8] Zaglmayr had applied 
the so-called scaled integrated Legendre polynomials in the construction, viz. 



n-l 



di, n>2, tG(0, 1], 



(29) 



where In (x) is the n-th order Legender polynomial. 



4.1 Edge functions 

For the completeness of our basis construction, in this subsection we record the results in [8]. 
Associated with each edge the formulas for these functions are given as 

'^Xm] = ^^^^^ ^ ^'^i - ^^^"^ ^ ^30) 

for the lowest-order approximation and 

= (^j+2(7e,,Afc2+AfcJ), j = 0,---,p-l (31) 
for higher-order approximations. 



4.2 Interior functions 

The interior functions are further classified into two categories: edge-based and bubble interior 
functions. By construction the normal component of each interior function vanishes on either edge 
of the reference 2-simplex K^, viz. 



n-..$t = 0, i = {1,2,3}, (32) 
where n^^ is the unit outward normal vector to edge ej. 

Edge-based interior functions: 



The tangential component of each edge-based function does not vanish on the associated only 
edge efc := [ki,k2] but vanishes the other two edges, viz. 

^'^^ • = 0' (33) 
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where t^^ is the directed tangent along the edge ej := [ji,j2]- The fohowing basis functions are 
proposed here as 



t,i 

e[fci,fc2] 



(0,2) 



2Ak 



k2 



1 - Afc, 



(34) 



where the function pj^^''^\») is the classical un-normalized Jacobi polynomial of degree i with a 
single variable [10], and the scaling coefficient is given by 



Ci = V2(i + 2)(i + 3)(2i + 3)(2i + 5), i = 0, 1, • • • ,p - 2. (35) 
The following orthonormal property of edge-based interior functions can be proved 

< ^e'E-fe]' ^e[fci,fc2] > = ^rnn, {m, n} = 0, 1, • • • , p - 2. (36) 

Interior bubble functions : 

The interior bubble functions vanish on the entire boundary dK^ of the reference 2-simplex K^. 
The formulas of these functions are given as 



= c^,nAoAiA2(i - Ao)™p42'2) (t::^) ^i'™""'"'^ (2^0 - 1) i^, i = 1,2, 



(37) 



where 



Crt 



(m + 3)(m + 4)(2m + 5)(2m + n + 6)(2m + n + 7)(2m + 2ra + 8) 



(m + l)(m + 2)(n + l)(n + 2) 



and 



< {m, n}, m + n < p — 3. 
One can again prove the orthonormal property of the interior bubble functions 



(38) 



where 



< {mi,m2,ni,n2}, mi + ni,m2 + ^2 < p - 3, {i, j} = 1,2. 

The following table shows the decomposition of the space (Pn(/^))^ for the ?{(div)-conforming 
triangular element. 



Decomposition 


Dimension 


Edge functions 


3(n+l) 


Edge-based interior functions 


3{n - 1) 


Interior bubble functions 


(n-2)(n- 1) 


Total 


(n + l)(n + 2) 



Table 1: Decomposition of the space (P„(X)) for the 'H(div)-conforming triangular element. 
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5 Basis functions for the tetrahedral element 



Our constructions are motivated by the work on the construction of ?^(div)-conforming hierarchi- 
cal bases for tetrahedral elements [7j. We construct shape functions for the ?^(div)-conforming 
tetrahedral element on the canonical reference 3-simplex. The shape functions are grouped into 
several categories based upon their geometrical entities on the reference 3-simplex [7]- The basis 
functions in each category are constructed so that they are orthonormal on the reference element. 

Any point in the 3-simplex is uniquely located in terms of the local coordinate system 
(C) ^) 0- The vertexes are numbered as vo(0, 0, 0), vi (1, 0, 0), V2(0, 1, 0), V3(0, 0, 1). The barycentric 
coordinates are given as 

Ao := 1 - e - - C, Ai := e, As := r?, A3 := C- (39) 
The directed tangent on a generic edge = [ji , js] is defined as 

.^^[hM (40) 

The edge is parameterized as 

7e, := Aj2 - Aj,, ji<i2. (41) 
A generic edge can be uniquely identified with 

ej := [ii, J2], Ji = 0, 1, 2, ji < j2 < 3, j = ji + 32 + sign(ii), (42) 

where sign(O) = 0. Each face on the 3-simplex can be identified by the associated three vertexes, 
and is uniquely defined as 

fji := b2,j3,i4], < {ii,j2,i3, ^4} < 3, j2 < is < J4- (43) 

The standard bases in M" are noted as Cj, z = 1, • • • , n, and n = {2, 3}. 
5.1 Face functions 

The face functions are further grouped into two categories: edge-based face functions and face 
bubble functions. 

Edge-based face functions: 

These functions are associated with the three edges of a certain face fj^, and by construction 
all have non-zero normal components only on the associated face fj^ , viz. 

n*^- ^SC..] = 0' ^^^Ji' (44) 

where n^-"= is the unit outward normal vector to face fj^.. 

The edge-based face functions for higher order have been proposed in [7] as follows: 

= ^i(7eJAfc, VAfc, X VAfc3, i = 0, • • • ,p - 1. (45) 
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For instance, for the face opposite to the vertex vo(0, 0, 0), fo := [1, 2, 3], the face functions related 
to edge e[l, 2] are given by 



^e[l,2] = ^i(^3 - A2)AiVA2 X VA3, i = 0, • • • ,p - 1. (46) 

However, it can be easily checked that the basis function given in (j45p in fact are not independent 
for p = 2 (for example, the sum of all the basis on a given face in fact equals to zero) and thus the 
proposed basis function is not complete. To remedy this degeneracy, two types of constructions of 
new hierarchical high-order independent edge-based face functions will be presented here. 

• Type One: high-order independent edge-based face functions 

In [9], the following orthonormal basis functions are given as 

- - A,)-^'-' ( - 1) ^1^. (47, 

Ci = V3(2i + 4)(2i + 5), i = 0,l,--- 



where 
and 



^"1 = {32, js}, h = {JsJa}, h < k2, ks = {j2, J3,i4} \ {ki, k2}. 
One can prove the orthonormal property of these edge-based face functions 

< ^e[Cfc2] ' >\K^=^mn, {m, n} = 0, 1, • • • , p - 1, (48) 

where 5mn is the Kronecker delta. Note that with our new construction, the edge-based face 
functions are all linearly independent, which is also verified by the fact that in the spectrum of the 
mass (Gram) matrix, none of the eigenvalues is zero. 



• Type Two: high-order independent edge-based face functions 



An alternative approach using the idea of recursion from [7] can also be used to construct 
independent edge-based face functions as follows. 

For p = 1, for each edge we have one face function for this edge as proposed in [7] 

5e[V°fc.] = ^^^iVA,,xVA,3, (49) 
and for p = 2, one additional new basis function can be constructed as 

%'fe] = ^^^i^^^2VAfc3xVAfc,, (50) 
which can be shown to satisfy the condition (I44p . and and for p > 3, the basis functions are given 

by 

- ^«(7eJ$|;;.,] + ^.-l(7ej5>5V°,,] (51) 

= ^i(7ej [AfciAfcjVAfcg X VAfcJ + ^i_i(7ej [Xk^VXk^ x VA^g] , i = 1, • • • ,p - 2. 
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It can be shown again numerically that there are exactly p functions that are independent and only 
whose normal component is non-zero only on the associated edge e^. 



Face bubble functions : 

The face bubble functions which belong to each specific group are associated with a particular 
face fj^. They vanish on all edges of the reference 3-simplex K^, and the normal components of 
which vanish on other three faces, viz. 

nf^fe.$^i„ = 0, jfc/ii. (52) 

The explicit formula is given as 



= .(1-A,,)-(1-A,,-A,3)"P(?"+^'^) - l) ( 



,j2 / \1 - Aj2 - Aj3 J |VAj3 X VAj^l ' 

(53) 
where 

i = C^Xj^Xj^Xj^, (54) 

where 

^„ _ yj{2n + 3)(m + n + 3)(m + 2n + 4)(m + 2n + 5)(2m + 2n + 7)(2m + 2n + 8)(2m + 2n + 9) 

"~ V(m + l)(m + 2y ' 

(55) 

and 

< {m, n}, m + n < p — 3. (56) 

By construction the face bubble functions share again the orthonormal property on the reference 
3-simplex K^: 

f f 

< ^>m\,ni, ^>m2,ri2 > \k^ = (^mima^nma, < {mi, 7712, TT-i, 712}, mi Tli , m2 + 71-2 < p - 3. (57) 

5.2 Interior functions 

The interior functions are classified into three categories: edge-based, face-based and bubble interior 
functions. By construction the normal component of each interior function vanishes on either face 
of the reference 3-simplex , viz. 

nf^.$t = 0, j = {0,l,2,3}. (58) 



Edge-based interior functions: 

The tangential component of each edge-based function does not vanish on the associated only 
edge efc := [ki, but vanishes all other five edges, viz. 

where t'^^ is the directed tangent along the edge := [ji,j2]- The shape functions are given as 
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r 



(60) 



where 



r , ; (2i + 4)(2i + 5)(2i + 7) 

Cj = (z + 3)^/ — — , ^ = 0, 1, • • • ,p - 2. 

? + 1 



Again one can prove the orthonormal property of edge-based interior functions: 

< ^eE.fca]' ^e[fci,fc2] > Ik'^ = ^mn, {m, n} = 0, 1, • • • , p - 2. (61) 

Face-based interior functions : 

These functions which are associated with a particular face ij^ have non-zero tangential com- 
ponents on their associated face only, and have no contribution to the tangential components on 
all other three faces, viz. 

n^^'^ X <^l^% = 0, jk^j^. (62) 
Further each face-based interior function vanishes on all the edges of the 3-simplex K^, viz. 

r^^- • f *S = 0. (63) 
The formulas of these functions are given as 



T 



(64) 

where l is given in (j54p and < {m,n},m + n < p — 3. The face-based interior functions enjoy the 
orthonormal property on the reference 3-simplex K^: 

< ^muni,^m2!n2 > Ik^ = ^mima'^nma,* = {1,2},0 < {iTll , m2 , ni , n2}, mi + ni,?n2 + 77,2 < p - 3. 

(65) 

Interior bubble functions : 

The interior bubble functions vanish on the entire boundary dK"^ of the reference 3-simplex K^. 
The formulas of these functions are given as 



^i;,. = ^^"^^'^^ (2Ai - 1) - l) [ ,_f;_^^ - l) ^ = 1, 2, 3, 

(66) 
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where 
where 
where 



X = C'€,m,nAoAlA2A3(l — Ai)™(l - Ai — A2)^ 



c] 



' + 2m + 2n + 9)(£ + 2m + 2n + 10)(2£ + 2m + 2n + ll)(m + 2n + 6) 



+ l)(m+l)(n+l) 



^2 



'(m + 2n + 7)(2m + 2n + 8)(n + 3)(n + 4)(2n + 5) 



+ 2)(m + 2)(n + 2) 



and 



< {I, m, n}, ^ + ?n + n<p — 4. 
Again, one can show the orthonormal property of the interior bubble functions 



where 



L c t e ■ 



< {£l,^2,"^l,"^2,?^-l,"-2},^l + mi + ni, ^2 + "^2 + "-2 < P - 4, {i, j} = 1, 2,3. 



In the following table we summarize the decomposition of the space {^n{K)f for the n{dW)- 
conforming tetrahedral element. 



Decomposition 


Dimension 


Edge-based face functions 


12n 


Face bubble functions 


2(n - 2){n - 1) 


Edge-based interior functions 


6(n - 1) 


Face-based interior functions 


4(n - 2){n - 1) 


Interior bubble functions 


(n-3)(n-2)(n- l)/2 


Total 


(n-M)(n + 2)(n + 3)/2 



Table 2: Decomposition of the space (P„(-ftr))^ for the 'H(div)-conforming tetrahedral element. 

6 Global divergence-free approximation 

For a p-th order polynomial approximation Up G ?^(div) and to ensure its global divergence- 
free condition, the idea is to add an higher-order interior bubble function Xb as an extra basis 
function, say 1 = max(m,p -|- 1) for the simplexes (m = 3 for triangular element and m = 4 for 
tetrahedral element) and q = max(2,p -|- 1) for rectangular and hexahedral elements. Namely, the 
new approximation in the augmented approximation space is 



ul = up + C^xl 



(67) 
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Due to the vanishing normal components of the interior bubble function xl along element bound- 
aries, the modified approximation u* still belongs to the space H^dW). 
Meanwhile, to satisfy the global divergence-free condition 

V • B = (68) 

for the Maxwell equation or 

V • u = (69) 
for the incompressible fluid flow, one can simply impose the condition 

V-(np + C^4) = 0, (70) 

which will give the unique value for the unknown coefficient C^. 

Now the piecewise divergence- free magnetic field u* € T-L{dW) will have zero spatial magnetic 
charge density (which equals to the divergence of the magnetic field by the Gauss' law for the 
magnetic field) and also zero surface magnetic charge density (which equals to the jump of the 
normal component of the magnetic field by the boundary condition of the magnetic field on material 
interfaces) along element interfaces, thus, artificial magnetic charge effect is avoided in u*. 

Finally, for the triangular element, the interior bubble function for degree p is given in equation 
(j37|) . For the tetrahedral element, the interior bubble function of degree p is given in equation (|66p . 
For the quadrilateral element, one can construct the interior bubble function of degree p as 

= Lp+2(2C - l)L2(2r/ - 1)1^ + ^2(2^ - l)Lp+2(2r? - l)j^. (71) 

For the hexahedral element, one can construct the interior bubble function of degree p as 

X?^^ = Lp+2{Cl)L2{V2)L2{C3)k + L2{Cl)Lp+2{V2)L2{C3)]r, + L2 (6 (??2)Lp+2 (C3)kc ■ (72) 



7 Conditioning of matrices 

The purpose of this section is twofold. Firstly, we check numerically that the newly constructed 
basis functions for 'H(div)-conforming triangular and tetrahedral elements are linearly independent, 
which is manifested by the fact that for each particular approximation order up to degree four, the 
condition number of the corresponding mass matrix is finite. Secondly, we want to show that for 
the approximation up to order three both the mass and stiffness matrices are reasonably well- 
conditioned. 

The components of the mass matrix are defined as 

M^,,^, :=< > l^rf, d = 2,3. (73) 

The mass matrix M is symmetric and positive deflnite, and therefore has real positive eigenvalues. 
The condition number of a real symmetric positive definite matrix A is calculated by the formula 

<A) = (74) 

where Amax and Amin are the maximum and minimum eigenvalues of the matrix A, respectively. 
For the incompressible fluid flows, e.g., governed by the Navier-Stokes equations [12] or by the 
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magnetohydrodynamics equations [13] the authors \12\ [T3] have appUed the mixed finite element 
for the spatial discretization. In particular, they \12\ I13j have used the ?^(div)-conforming element 
for the Laplacian Au of the velocity u. In this case, we have the stiffness matrix S, which is defined 
component- wise as 

Se„e, := < V^e, : V$,, > d = 2, 3. (75) 

The stiffness matrix S is symmetric and semi-positive definite, and therefore has real non-negative 
eigenvalues. The condition number of the stiffness matrix S is calculated by the formula (j74p with 
the zero eigenvalue excluded. 

With the triangular element and for the polynomial approximations p = {1,2,3,4}, the condi- 
tioning is summarized in Table [3j From the table we can see that the condition number is bounded 



Order p 


Mass 


Stiffness 


1 


2.016el 


1.040el 


2 


8.804el 


5.959el 


3 


9.847e2 


4.197e2 


4 


1.286e4 


8.843e3 



Table 3: Condition numbers of the mass matrix M and stiffness matrix S from the basis for the 
■H(div)-conforming triangular element. 

for each order of approximation. Moreover, up to the fourth order, the mass and stiffness matrices 
are both well conditioned. 

With the tetrahedral element and for the polynomial approximations p = {1,2,3,4}, the con- 
dition numbers of the mass matrix are shown in Table HI 



Order 
P 


Mass 


Stiffness 


Ratio 


Type 1 


Type 2 


Type 1 


Type 2 


Mass 


Stiff. 


1 


3.084el 


3.084el 


1.989el 


1.989el 


l.OOOeO 


l.OOOeO 


2 


6.987e3 


7.733e4 


3.395e3 


5.917e4 


0.090e0 


0.057e0 


3 


3.412e6 


2.138e6 


1.094e6 


9.919e5 


1.491e0 


0.919e0 


4 


5.972e9 


2.639e7 


2.883e9 


2.289e7 


2.198e2 


1.215e2 



Table 4: Condition numbers of the mass matrix M and stiffness matrix S from the bases with two 
different types of edge-based face functions for the ?^(div)-conforming tetrahedral element. 

Again from this table we see that the condition number is finite for each order of approximation. 
Further up to the third order, both the mass and stiffness matrices are well conditioned. For order 
p = 2 the conditioning is better with the Type-one basis while for p = 4, the conditioning is better 
with the Type-two basis. For the third order p = 3, the performance with both types of bases is 
about the same. 
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8 Concluding remarks 



In this paper we focus our attention on hierarchical ?^(div) basis functions for solving the magneto- 
hydrodynamics (MHD) equations numericahy so that the divergence-free condition on the magnetic 
field is rigorously guaranteed. The idea is to use an interior bubble function from the proposed high- 
order hierarchical basis as the additional freedom to impose the divergence-free Gauge condition 
for the magnetic field. We have summarized four bases for the 7^(div)-conforming elements, viz. 
the quadrilateral and triangular elements for 2-D and the hexahedral and tetrahedral elements for 
3-D. The linear independence of the basis functions for the two simplicial elements has numerically 
been checked. Good matrix (mass and stiffness) conditioning has also been shown up to the fourth 
order for 2-D and up to the third order for 3-D. Further work will include the implementation of 
the proposed divergence-free basis to solve the magnetohydrodynamics (MHD) equations in 2-D 
and 3-D. 
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